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We investigate the nonlinear wave dynamics of origami-based metamaterials composed of Tachi-Miura poly¬ 
hedron (TMP) unit cells. These cells exhibit strain softening behavior under compression, which can be tuned 
by modifying their geometrical configurations or initial folded conditions. We assemble these TMP cells into a 
cluster of origami-based metamaterials, and we theoretically model and numerically analyze their wave trans¬ 
mission mechanism under external impact. Numerical simulations show that origami-based metamaterials can 
provide a prototypical platform for the formation of nonlinear coherent structures in the form of rarefaction 
waves, which feature a tensile wavefront upon the application of compression to the system. We also demon¬ 
strate the existence of numerically exact traveling rarefaction waves. Origami-based metamaterials can be highly 
useful for mitigating shock waves, potentially enabling a wide variety of engineering applications. 

PACS numbers: 45.70.-n 05.45.-a 46.40.Cd 


I. INTRODUCTION 

Recently, origami has attracted a significant amount of at¬ 
tention from researchers due to its unique mechanical prop¬ 
erties. The most evident one is its compactness and deploya¬ 
bility, which enables various types of expandable engineering 
structures, e.g., space solar sails GLU and solar arrays H. Bi¬ 
ological systems also exploit such compact origami patterns, 
such as foldable tree leaves for metabolic purposes 0] and 
stent grafts 0. Another useful aspect of origami-based struc¬ 
tures is that origami patterns can enhance static mechanical 
properties of structures. For instance, structural bending rigid¬ 
ity for thin-walled cylindrical structures can be significantly 
improved by imposing origami-patterns 0. These origami 
patterns are used not only for space structures, but also in 
commercial products (e.g., beverage cans) in order to reduce 
the thickness of thin-walled structures without sacrificing their 
buckling strength 0. 

Within the considerable progress made in the mechanics of 
origami-based structures, however, the primary focus has been 
placed on the static or quasi-static properties of origami. For 
example, recent studies attempted to fabricate origami-based 
metamaterials with an eye towards investigating the deploy¬ 
able, auxetic, and bistable nature of origami structures ]8|- 
m. Fimited work has been reported on the impact response 
of origami-based structures [12], and their wave dynamics 
is relatively unexplored. Plausibly, this lack of studies on 
the dynamics of origami-based structures can be attributed 
to the intrinsic characteristic of typical origami structures, 
which exhibit limited degrees of freedom (DOF) during their 
folding/unfolding motions. This is particularly true for rigid 
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origami, in which the deformation takes place only along 
crease lines while origami facets remain rigid in dynamic con¬ 
ditions. The rigid origami features single-DOF motions ide¬ 
ally, and thus, the studies on their wave dynamics have been 
more or less absent under this rigid foldability assumption. 

In this study, we use a single-DOF rigid origami structure 
as a building block to assemble multi-DOF mechanical meta¬ 
materials, and analyze their nonlinear wave dynamics through 
analytical and numerical approaches. Specifically, we employ 
the Tachi-Miura polyhedron (TMP) lfl3l ITill as a unit cell of 
the metamaterial as shown in Fig. [Q The TMP cell is made 
of two adjoined sheets (Fig. |T} a)), and changes its shape from 
a vertically standing planar body to a horizontally flattened 
one while taking up a finite volume between the two phases 
(Fig. Eh))- This volumetric behavior is in contrast to conven¬ 
tional origami-patterns that feature planar architectures and 
in-plane motions (e.g., Miura-ori sheets fl5tl ). In this study, 
we first characterize the kinematics of the TMP cell, show¬ 
ing that it exhibits controllable strain-softening behavior. By 
cross-linking these TMP unit cells in a horizontal layer and 
stacking them up vertically with separators, we form a multi- 
DOF metamaterial as shown in Fig. |TJc). We then conduct an¬ 
alytical and numerical studies to verify that these multi-DOF 
origami structures can support a nonlinear stress wave in the 
form of a so-called rarefaction wave, owing to the strain soft¬ 
ening nature of the assembled structure. 

The rarefaction wave, which can be viewed as an acoustic 
variant of a depression wave 0, has been studied in various 
settings, including systems of conservation laws | fl7tl . More 
recently, it was proposed in the_context of discrete systems 
with strain-softening behavior ifUil fl9tl . Interestingly, these 
rarefaction waves feature tensile wavefronts despite the ap¬ 
plication of compressive stresses upon external impact (see 
the conceptual illustrations in Fig. |T|c)). In that light, they 
are fundamentally different from the commonly encountered 
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FIG. 1: (Color online) (a) Flat front and rear sheets of the TMP 
with mountain (solid lines) and valley (dashed lines) crease lines, (b) 
Folding motion of the TMP unit cell, (c) System consisting of TMP- 
based metamaterials and rigid separators stacked vertically. Each 
layer consists of nine inter-linked TMP unit cells. Conceptual illus¬ 
trations of incident compressive waves and transmitted rarefaction 
waves are also shown. 


dynamical response of nonlinear elastic chains which sup¬ 
port weakly or even strongly nonlinear traveling compression 
waves I19 h 22I1 . Recently, in a quite different setting of tenseg- 
rity structures, rarefaction waves have been identified compu¬ 
tationally in the elastic softening regime 0. 

Our main scope within the present work is to verify the 
formation and propagation of rarefaction waves in origami- 
based metamaterials via two simplified models: a multi-bar 
linkage model and a lumped mass model. In both cases, we 
confirm that the origami structure disintegrates strong impact 
excitations by forming rarefaction waves, followed by other 
dispersive wave patterns to be discussed in more detail below. 
We also validate the nonlinear nature of the stress waves by 
calculating the variations of wave speed as a function of ex¬ 
ternal force amplitude. Notably, we observe the reduction of 
wave speed as the excitation amplitude increases, which is in 
sharp contrast to conventional nonlinear waves seen in nature 
or engineered systems 0. In the case of the lumped mass 
model, we find numerically exact traveling waves. We pro¬ 
vide a precise characterization of the wave speed and ampli¬ 
tude relationship and a way to evaluate the robustness of the 
rarefaction waves through dynamical stability computations. 
The findings in this study provide a foundation for building a 
new type of impact mitigating structure with tunable charac¬ 


teristics, which does not rely on material damping or plastic 
deformation. This study also offers a platform for exciting the 
rarefaction pulse - a far less explored type of traveling wave 
- and examining its characteristics in considerable detail. 

The Manuscript is structured as follows: In Sec. [TT] we de¬ 
scribe the two simple models of origami-based metamaterials: 
the multi-bar linkage model and the lumped mass model. In 
Sec. um we conduct numerical simulations of wave propaga¬ 
tion upon impact on the chain boundary and compare the wave 
dynamics obtained from these two models. Then, in Sec. IIVI 
we find numerically exact rarefaction waves of the lumped 
mass model. Lastly, concluding remarks and future work are 
given in Sec. El 


II. MODELING OF ORIGAMI-BASED STRUCTURES 
A. Multi-bar Linkage Model 

We begin by modeling a single TMP cell as shown in Fig. U 
For the sake of simplicity, we focus on the folding motion of 
two adjacent facets along the horizontal crease line as marked 
by the red line in Fig. EJa). Preserving the key features of the 
TMP, such as rigid foldability and single-DOF mobility, we 
can model the folding/unfolding motion of the origami facets 
into a simple ID linkage model as shown in Fig. EJb). Here, 
the unit cell consists of two rigid bars (mass m and length 
2 L), and the center-of-mass coordinates of those two bars are 
( z \, j/i, 6\) and (z 2 , yi, $ 2 )- The hinge that connects the two 
bars is equipped with a linear torsional spring with the torsion 
coefficient kg. The left end of the linkage structure is sup¬ 
ported by a roller joint, which is allowed to move only along 
the z-axis up on the application of external force F ex . The 
right end is fixed by a pin joint. Therefore, the inclined angle 
of the linkage, 6 1 , is the only parameter required to describe 
the motion of this unit-cell system. This corresponds to the 
single-DOF nature of the TMP cell. 

By using the principle of virtual power fi^l, we derive the 
following equation of motion (see Supplemental Material for 
details |25|): 

(mi 2 /2 + J/2 + 2mL 2 cos 2 0i) §1 — mL 2 d\ sin20i 

+k 9 (6>i - 0i, o ) = —F ex L cos 9 1 . (1) 

Here J is the bar’s moment of inertia (J = m ^~), and ()\.0 
is the initial folding angle (i.e., no torque applied at the hinge 
in this initial angle). In the quasi-static case (i.e., acceleration 
and velocity terms are much smaller compared to the exter¬ 
nal excitation and spring force terms), we obtain the force- 
displacement relationship as follows: 

pex = fee (#i - <9i,o) 

L cos 6\ 

Using Eq. (0 and the axial displacement expression u = 
4L(sin o — sin 0i), we can calculate the force-displacement 
response as shown in Fig. 0c). We observe that the system 
exhibits strain softening behavior in the compressive region. 
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FIG. 2: (Color online) (a) TMP unit cell, (b) Two-bar linkage model representing the folding motion of the two facets as marked in red 
lines in (a), (c) Force-displacement relationship of the TMP unit cell with L = 5 mm, kg = l.ONm/rad, and different initial folding angles: 
8i ,o = 30°, 45°, and 60°. Dashed line indicates a power law approximation of 6 * 1,0 = 4:5° case. 


whereas the system shows strain hardening response in the 
tensile domain. Also, it is interesting to find that this strain 
softening/hardening behavior can be tuned by controlling the 
initial folding angle, 0 1O . 


hi h 2 h N 



FIG. 3: Schematic illustrations of (a) Multi-bar linkage model and 
(b) Lumped mass model. 


Based on the kinematics of the single unit cell as ex¬ 
pressed in Eq. Q, we model a chain of (V-TMP cells as 
shown in Fig. [3a). In this model, each unit cell is con¬ 
nected by pin joints, which are allowed to move along the 
z-axis. Let the general coordinate of this A'-DOF system be 
q = [ 6*1 ••• 9j ••• 6 fv] . Given the identical initial angles 
imposed on the unit cells, the equation of motion for this sys¬ 
tem can be expressed as 


G T MGq + G T MGq = G T f ex (3) 


where 

M = diag [Mi • • • Mjv] 


Mi 

= diag [m m 

J m m J] , 
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Also, f ex is an external force vector defined as follows 

i T 


_ 


(ff 


(f rY 


(f W) 


(4) 


where 


[F ex , 0, - 2k e (6>1 - 01,o) - F ex L COS01, 

. 0 , 0 , —2kg ( 0 ljO — 0 i)] T if j = 1 

j 1 [0, 0, -2 ke (01 - 01,0), 

0, 0, —2ko (0j,o — 0j)] T if . 7 = 2 ...TV 

See Supplemental Material for the details of this deriva¬ 
tion [25]. 


B. Lumped Mass Model 

In this section, we introduce a lumped mass model, in 
which a chain of origami cells is modeled as lumped masses 
connected by nonlinear springs (see Fig. [3b)). The strain soft¬ 
ening behavior of the TMP unit cell considered herein leads 
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to the following power-law relationship: 

F ex = A5 n (5) 

where S is the compressive displacement, and the coefficient 
A and the exponent n are the constant values determined by 
curve fitting of Eq. ©. 

Since the power-law relationship in Eq. © assumes only 
a positive displacement as an argument, we need to apply a 
displacement offset (do) towards the tension side, so that the 
lumped mass model can approximate the force-displacement 
curve of the multi-bar linkage model not only in the com¬ 
pressive region, but also in the tensile domain. In Fig. ©c), 
the dashed curve shows the fitted power-law relationship for 
the multi-bar linkage model, where the black circle represents 
(along the horizontal axis) the displacement offset d 0 . 

By using this simple force-displacement relationship, we 
can derive a general expression of the equation of motion as 
follows: 


Muj — A [do + Sj-i,j\ + — A [do + Sjj+ 1 ]” (6) 

where M is the lumped mass corresponding to 2 to, n G R, 
and the [+] sign outside of the brackets indicates that we take 
only positive values of the strain Sjj+i = Uj — it ;+ i. Note 
that this form of equation has been used widely for analyzing 
nonlinear waves propagating in discrete systems in the case 
of strain-hardening interactions (i.e., n > 1 in Eq. ©, e.g., 
granular crystals). Therein, the formation and propagation of 
nonlinear wave structures, such as solitary waves fl9l 120(1 and 
discrete breathers ll2ll [22i l. have been well studied. The in¬ 
terpretation of origami dynamics via this nonlinear lumped 
mass system opens up a broad, novel potential vein of stud¬ 
ies. Indeed, one advantage of modeling the origami lattice in 
this way is that many tools and results obtained in the con¬ 
text of granular crystals can be applied in our setting. For ex¬ 
ample, the recent work of [ 18] examined a one-dimensional 
discrete system under the power-law relationship of strain¬ 
softening springs (i.e., n < 1 in Eq. ©). This study reported 
the propagation of rarefaction waves through dynamic sim¬ 
ulations and a long wavelength approximation, where it was 
shown that the width of the rarefaction wave is independent of 
the wave speed. Likewise, the analysis of nonlinear waves in 
post-buckled structures has been also attempted using a sim¬ 
ilar discrete system li26il . In this article, we extend the theo¬ 
retical results in such nonlinear-spring systems by introduc¬ 
ing a systematic tool for the computation of numerically exact 
traveling waves, which will be discussed in Sec.[IV] We also 
address the subject of their dynamical stability in the Supple¬ 
mental Material l25ll . 


m. NUMERICAL SIMULATIONS 

To examine the dynamic characteristics of the origami- 
based structure and compare the results from the two reduced 
models, we conduct numerical computations of wave propa¬ 
gation under a compressive impact. Also, we apply various 
amplitudes of impact force to the multi-bar linkage model in 


order to examine the speed of both compressive and tensile 
strain waves, especially focusing on the dominant traveling 
wave. 


A. Waveform analysis 


We perform numerical computations where a compressive 
impact is applied to the first unit cell with the right end of the 
7V-th unit cell kept fixed as shown in Fig. ©a). The strain 
waves propagating in a uniform chain of TV = 400 unit cells 
are examined numerically. In the case of the multi-bar linkage 
model, the relative strain is defined as 


m = 


foj ,0 hj 

h j, 0 


(7) 


where hj = ALsmdj and hj,o = AL sin 6jo (see Fig.©b)). 
The numerical constants used in the calculation are the follow¬ 
ing: L = 5mm, m = 0.39g, ke = l.ONm/rad, and 9j t o = 
45°. To apply impact excitation, we impose F ex = 100 N for 
the first 1 ms and F ex = 0 N after the first 1 ms in our sim¬ 
ulations. From the force-displacement curve based on these 
constants, we obtain n = 0.64 and A = 2,938N/m rl , given 
an initial displacement offset of do = 2.1 mm for the power- 
law approximation. In the case of the lumped mass model, the 
relative strain is defined as 


n, = u ’*\ ( 8 ) 

do 

Figures©a) and (b) show space-time contour plots of strain 
wave propagation under compressive impact, while Figs. [4] c) 
and (d) show the strain waveforms corresponding to t, 3,40, 
and 70 ms. After the impact force is applied to the system, the 
first compressive impact attenuates quickly as the strain waves 
propagate through the system, and then a rarefaction wave ap¬ 
pears in front of the first compressive wave (see the insets as 
well as the arrows (1) and (2) in Fig. ©>. It should be also 
noted that due to the strain-softening behavior, the amplitude 
of the compressive force is reduced drastically as the wave 
propagates along the chain. Since both the multi-bar linkage 
model and the lumped mass model have this strain-softening 
nature, the same type of rarefaction waves is observed. 

In addition, the inset of Fig. ©c) shows the magnified view 
of the leading edge of the propagating strain wave. This lead¬ 
ing wave is created due to the effect of inertia in the multi-bar 
linkage model. That is, when the first unit cell folds right 
after the compressive impact, the second unit cell is pulled 
by the first unit cell before the compressive force propagates 
to the next unit cell. Therefore, the tensile strain appears in 
front of the first compressive wave in the multi-bar linkage 
model. Comparing the numerical results of the two mod¬ 
els, the lumped mass model captures the multi-bar linkage 
model dynamics even quantitatively at short times, while the 
agreement between the two becomes qualitative at longer time 
scales. 

Let us also note in passing that in the wake of this primary 
rarefaction pulse, we observe radiative dispersive wavepack- 
ets both in the multi-bar linkage model and in the lumped mass 
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FIG. 4: (Color online) Space-time contour plots of strain wave propagation based on (a) the Multi-bar linkage model and (b) the Lumped mass 
model. Insets show the magnified view of rarefaction waves. Temporal plots of strain waves using (c) the Multi-bar linkage model and (d) 
the Lumped mass model. The inset in (c) shows the magnified view of the leading edge. The arrows (1) and (2) point to the rarefaction wave 
present in the dynamics. 


model. These wavepackets apparently travel maximally with 
the speed of sound in the medium, while the rarefaction pulse 
outrunning them is apparently supersonic. We will return to 
this point to corroborate it further by our numerical bifurcation 
analysis in the next section. Additionally, it should be noted 
that in the lumped mass model, highly localized structures 
with a clear envelope can be discerned (see e.g., the vicin¬ 
ity of unit number 150 of the 70 ms panel of Fig.|4fd)), which 
seem to have the form of breather excitations, which are ex¬ 
ponentially localized in space and periodic in time lf27l [28tl . 
A closer inspection of Fig. @[b) also seems to suggest that 
such coherent wavepackets travel more slowly than the dis¬ 
persive radiation. The multi-bar linkage model also exhibits 
such time-periodic patterns, but there is no clear signature of 
spatial localization. While these nonlinear wave structures are 
worth investigating, this topic is beyond the scope of this pa¬ 
per, and we do not explore them further here. 


B. Wave speed analysis 

The propagation speed of strain waves is now investigated 
numerically under various amplitudes of impact force. The 
wave speed is approximated as follows 


where ho is the initial height of the unit cell, and At is the 
time span in which the strain wave propagates from the first 
unit cell to the N-th unit cell (see Fig.[5[a)). The propagating 
wave speeds calculated are depicted in Fig. [5jb) under three 
different initial folding angles: 0 h o = 35°, 45°, and 55°. It 
is evident that the wave speed is altered by the impact force, 
which is one of characteristics of nonlinear waves. However, 
it should be noted that in the compressive regime, the wave 
speed decreases as the compressive impact increases. This is 
in sharp contrast to conventional nonlinear waves formed in 
the system of strain-hardening lattices EH. A different 
trend is observed in the tensile regime, where the wave speed 
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FIG. 5: (Color online) (a) Surface map of strain field to calculate wave speed, (b) Wave speed of strain waves as a function of external force 
ranging from —150 N to +150 N. Numerical simulations are based on L = 25 mm, rrij = 19.7 g, N = 20 and kg = 1.0 Nm/rad. 


increases as the tensile impact increases. It is also notewor¬ 
thy that the wave speed curve can be shifted by changing the 
initial folding angle. Therefore, we can control the speed of 
the waves propagating through the origami-based metamateri¬ 
als by altering their geometrical configurations, implying their 
inherent dynamical tunability. 

IV. EXACT RAREFACTION WAVES OF THE LUMPED 
MASS MODEL 

We now turn our attention to a more systematic analysis and 
understanding of the rarefaction waves in the simpler lumped 
mass model; notably, our conclusions here in that regard are of 
broader interest to previously discussed settings such as those 
of l HH l23l l. Based on the previous analysis, we numerically 


investigate the existence and dynamical stability of exact rar¬ 
efaction waves of the lumped mass model [cf. Eq. ©]. In 
particular, we consider the model in the strain variable 5jj +1 
written as 

1 = A{[d 0 + — 2[do + 5j,3+i]+ 

+ [do + ^j+ij+ 2 ]” }• ( 10 ) 

The existence and the spectral stability of traveling waves of 
Eq. ( ITOl ) with wave speed c must be examined through the 
ansatz 6jj + i(t) = S(j — cf) := $(£,f), i.e., going to the co¬ 
traveling wave frame where the relevant solution appears to 
be steady and hence amenable to a spectral stability analysis. 
Then, $ solves the advance-delay differential equation 


f) — —c * 2 <E»^(^, f) + 2c<b{ t (£, f) + — | [do + $(£ — 1, i)]+ ^ 2 [do + 3>(£, f)] + + [do + 3>(£ + 1, f)] + j- (11) 


Traveling waves of Eq. (fToh correspond to stationary (time 
independent) solutions <3?(£, t) = 0 (f) of Eq. i l 1 I | i. satisfying 

0 = -c 2 ^ + ^{[do + ^-l)]”-2[d o + 0(O]" 

+ [do + </>(£ + !)]+}• (12) 

To obtain numerical solutions of Eq. (ITH i. we employ a uni¬ 
form spatial discretization of £ consisting of l points £*, (k = 
— 1=2,..., 0,..., 1=2) with lattice spacing A£ chosen such 
that q = 1/A£ is an integer. Then, the field 0(£) is replaced 
by its discrete counterpart, i.e., </+ := </>(£*,) = </>(fcA£). 
The second-order spatial derivative appearing in Eq. ([Til l 
is replaced by a modified central difference approximation 
(4>k- 2 — 2 4>k + <f>k- |_2)/(4A£ 2 ). The reason for this choice of 


central difference is connected to the stability calculation to be 
discussed in the Supplemental Material 1251 1. Using this dis¬ 
cretization, Eq. dTH) becomes the following root-finding prob¬ 
lem, 


0 = 


2 — 2 <t>k + <t>k+2 


4A £ 2 

— 2 [dg + 4>k\+ + [do + </>fc+qr]_|_} 




<\>k- 


q\+ 


(13) 


which is solved via Newton iterations. We employ periodic 
boundary conditions at the edges of the spatial grid. We are in¬ 
terested specifically in rarefaction waves, and thus we use the 
profiles obtained via the numerical simulations of Sec. IIII Al 
to initialize the Newton solver, see e.g. arrow (2) of Fig. [4j d). 
Herein, we consider an origami lattice with L = 25 mm, 
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kg = l.ONm/rad and 9 = 55°. The corresponding best- 
fit values of the parameters of the lumped-mass model are 
A = 280N/m n , n = 0.53, m = 19.7g with M = 2m and 
do = 12 mm. 

In Fig. [6j numerically exact rarefaction waves (i.e., solu¬ 
tions of Eq. ( IT3l > with a prescribed tolerance) are presented 
for various values of the wave speed c. In particular, Fig. |6(a)| 
shows the rarefaction waves in terms of the relative strain vari¬ 
able <fr/do, while Fig. |6(b)| shows the corresponding relative 
momenta (j)'/do- Note that the tails decay to zero monoton- 
ically, implying that the traveling structure does not resonate 
with the linear modes of the system, as the wave is super¬ 
sonic. It is not surprising then that our parametric continu¬ 
ation in the wave speed c reveals a critical minimum value 

c s = nAd/j^ 1 /M = 173.5 m/s, which is the sound speed 
of the chain (see the vertical dashed-dot gray line of Fig. |6(c)l i. 
This is consistent with the long-wavelength analysis of ll 1811 
and also with our observations of the previous section indicat¬ 
ing that the wave outruns the small amplitude radiation tails 
behind it. Thus, similarly to systems with n > 1 EEi, 
the rarefaction waves of the origami lattice are traveling faster 
than any linear waves of the system. However, in contrast 
to solitary waves in systems with n > 1, the amplitude of 
the rarefaction waves in the origami system have a natural 
bound determined by the precompression do of the system, 
in which case the particles come out of contact (see the hor¬ 
izontal dashed black line of Fig. |6(c)| >. Although waves with 
amplitude exceeding this value are in principle possible, we 
were unable to identify any ones such numerically. An inter¬ 
esting open problem would be to prove rigorously if such a 
bound exists. Another interesting related problem is if there 
is a critical maximum value of c. Our numerical continuation 
algorithm did indeed terminate due to lack of convergence at 
c « 201.6 m/s, but this could have been a result of the ill- 
conditioned nature of the Jacobian matrix as the amplitude 
approached the critical limit of do. 

The robustness of a solution <fp of Eq. (12) can be inves¬ 
tigated through a spectral stability analysis. To that end, we 
substitute the linearization ansatz $(/;,/) = <p° + ea(£)e At 
into Eq. (11), which yields an eigenvalue problem at order 
e (see the Supplemental Material |25|). We considered so¬ 
lutions at various wave speeds c and found in each case at 
least one eigenvalue with a small real part, indicating a (very 
weak) instability. However, the eigenvalues are highly sensi¬ 
tive to e.g. lattice size and choice of discretization, suggesting 
that these instabilities may be “spurious”. To check this, we 
performed dynamical simulations of the perturbed rarefaction 
waves at the level of Eq. ( ITOt and found that they are all robust 
against small perturbations (see Fig. 8 of the Supplemental 
Material lf25h ). This suggests that the very weak instabilities 
predicted via the spectral stability analysis are indeed spuri¬ 
ous. While a heuristic argument for the presence of spurious 
instabilities is provided in the Supplemental Material B25I1 . the 
construction of a mathematically consistent algorithm for the 
computation of eigenvalues in this context remains an impor¬ 
tant open problem. 


V. CONCLUSIONS & FUTURE CHALLENGES 

In the present work, we investigated nonlinear wave dy¬ 
namics in origami-based metamaterials consisting of building 
blocks based on Tachi-Miura polyhedron (TMP) cells. We 
analyzed the kinematics of the TMP unit cell using a sim¬ 
ple multi-bar linkage model and found that it exhibits tun¬ 
able strain-softening behavior under compression due to its 
geometric nonlinearity. We observed that upon impact, this 
origami-based structure supports the formation and propaga¬ 
tion of rarefaction waves. The resulting evolution features a 
tensile wavefront despite the application of compressive im¬ 
pact. A further reduction was also offered based on the fitted 
force-displacement formula for a single cell, in the form of 
a lumped mass model. In the latter case we obtained numeri¬ 
cally exact rarefaction waves and studied their spectral and es¬ 
pecially dynamical stability. The dynamical features observed 
herein may constitute a highly useful feature towards the ef¬ 
ficient mitigation of impact by converting compressive waves 
into rarefaction waves and disintegrating high-amplitude im¬ 
pulses into small-amplitude oscillatory wave patterns. We 
also demonstrated the potential tunability of the wave speed 
by altering initial folding conditions of the origami-based 
structure, which naturally opens up the feasibility of control¬ 
ling stress wave propagation in an efficient manner. 

The rather unique nonlinear wave dynamics of origami 
structures can lead to a wide range of applications, such as 
tunable wave transmission channels and deployable impact 
mitigating layers for space and other engineering applications. 
On the theoretical/computational side, there is also a large 
number of intriguing questions that are emerging. For one, 
a more detailed comparison of the coherent structure propa¬ 
gation in the multi-bar linkage model vs. that of the lumped- 
mass model would be an interesting topic for further consider¬ 
ation. This would help uncover the dynamical features leading 
to the apparent weak amplitude decay in the former, while the 
latter contains robust solutions and sustained long-time propa¬ 
gation. Still at the single wave level, an exploration of the del¬ 
icate issues of spectral stability by means of different numeri¬ 
cal methods and of the corresponding dynamical implications 
would be of particular interest. Subsequently, understanding 
further the dynamics and interactions of multiple rarefaction 
wave patterns would also be a relevant theme for future inves¬ 
tigations. These topics are currently under active considera¬ 
tion and will be reported in future publications. 
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FIG. 6: (Color online) Summary of numerical results on continuations of rarefaction waves over wave speed c with l = 4001 points and 
= 1/13: (a) Relative strain profiles for various values of the wave speed c. (b) Relative momenta corresponding to (a), (c) Maximum of 
the absolute value of the relative strain variable as a function of the wave speed. Note that the horizontal dashed black line corresponds to the 
value of pre-compression in normalized units (or, equivalently, do in physical units), while the vertical dashed-dot gray line corresponds to the 
value of the speed of sound c a of the medium. 
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